Chandrasekhar Theory of Ellipsoidal Electromagnetic Scatterers 
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A number of new problems in remote sensing and identification of buried compact metallic targets 
motivate the search for new models that, if not exact, at least enable extremely rapid numerical 
predictions of electromagnetic scattering/induction data. Here the elegant Chandrasekhar theory 
of the electrostatics of charged ellipsoids is used to develop an essentially exact, extremely efficient 
description of low- to intermediate frequency (or late- to intermediate-time) responses of ellipsoidal 
targets. Comparisons with experimental data demonstrate that, together with a previously de- 
veloped theory of the high frequency (or early time) regime, the results serve to cover the entire 
dynamic range encountered in typical measurements. 



Exactly soluble models in the theory of EM propa- 
gation and scattering are essentially limited to horizon- 
tally stratified or spherically symmetric geometries, with 
limited results also available for cylindrically symmetric 
waveguide geometries However, there are a num- 
ber of new problems in remote sensing and classification 
of buried compact metallic targets that require a wider 
class of solutions that, if not exact, at least support rapid 
numerical evaluation. These include a number of long- 
standing economic and humanitarian problems, such as 
clearance of unexploded ordnance (UXO) from old test 
ranges. The most difficult technological issue is not the 
detection of such targets, but rather the ability to dis- 
tinguish between them and harmless clutter items, such 
as pieces of exploded ordnance. Since clutter tends to 
exist at much higher density, even modest discrimination 
ability leads to huge reductions in remediation costs. 

The most promising technologies, least influenced by 
the complexities of the heterogeneous earth, are essen- 
tially low frequency (< 10 kHz) metal detectors that emit 
a series of sharply terminated, well spaced pulses, and 
measure the decaying induced currents during the (~ 25 
ms) quiet interval between pulses. At these low frequen- 
cies, detailed spatial resolution is lost [ij, but the induced 
voltage curve V{t) has substantial structure, and one 
may hope to extract target geometry information from 
it. This inference, however, is indirect fsl and solution of 
this inverse problem requires careful comparison of pre- 
dicted and measured responses. This requires accurate 
models that go well beyond existing exact solutions. 

A number of classes of buried objects, especially UXO, 
have a cylindrical geometry. As a step in the direction 
of modeling such objects, I consider here EM scatter- 
ing from ellipsoids, with arbitrary axes a = (01,09,03). 
Although fully analytic solutions are not possible [j], it 
is demonstrated here that Chandrasekhar's elegant ap- 
proach to the electrostatics of heterogeneously charged 
ellipsoids 0] enables the major part of the computation 
to be performed analytically, reducing it to a matrix di- 
agonalization problem. The size of the matrix grows as 
one pushes earlier in time, closer to pulse termination 
(where the high frequency components of the pulse spec- 



trum have not yet had a chance to decay away), but en- 
ables an essentially exact description at intermediate- to 
late-time at minimal numerical cost. As will be demon- 
strated through comparisons with laboratory data on ar- 
tificial spheroidal targets, combining this approach with a 
complementary early-time approach Q enables accurate 
predictions at all time (or frequency) scales. 

For clarity of exposition, the theory will be presented 
for nonmagnetic target and background, fj, = fib, where 
fib is the (uniform) background permeability. Magnetic 
targets exhibit both a conductivity and a permeability 
contrast, which complicates the analysis but does not 
affect the basic approach. Details of this generalization 
will be presented elsewhere. 

For nonmagnetic systems, the Maxwell equations may 
be reduced, in the frequency domain, to a single equation 
for the electric field 



V X V X E - K^E = S, 



(1) 



where = efik^, k = lj/c, and source term S ~ 
{A'Kifik/cy^s where j 5 is the source current density, rep- 
resenting in this case the transmitter coil. Inside the 
metallic target, at all frequencies of interest here, the 
dielectric function is dominated by the dc conductiv- 
ity, e ~ Airia/uj. The background may be insulating or 
weakly conducting, but we work in the very high contrast 
limit le/ebl >> 1 [7]. 

To obtain a closed equation, restricted to the finite 
target domain Vc, we use the Green function approach. 
A background field Ef, is defined by solving ^ for the 
same source S, but in the absence of a target, = k^. 
By subtracting this equation from ([T]), and applying the 
background tensor Green function G, defined by 

V X V X G(x, x') - KbG(x, x') = l(5(x - x'), (2) 

where 1 is the 3x3 identity matrix and the curls act on 
the first index of G, (H]) takes the integral form 



E(x) - Eb(x) 



d3x'Q(x')G(x,x') •E(x'), (3) 



where the contrast function Q — — vanishes out- 
side the target volume Vc- For a uniform background 



2 



one obtains G = (1 + K^'^W)g, where g(x, x') = 
giKb|x-x I ^47r|x — x'l is the scalar Helmholtz Green func- 
tion. Restricting x € 14, (|3]) becomes an equation for the 
internal field Eint alone. Given Ei„t(x), the external field 
Eext(x), X ^ T4, follows by direct integration. 

Over the measurement domain outside Vc, Hb can vary 
strongly, encompassing soil, rock, surface plants, air, etc. 
The full E field is therefore unpredictable. However, an 
EM induction (EMI) measurement is sensitive only to 
the "magnetic" contribution to E. Specifically, one may 
divide E into curl-free and divergence free contributions, 
and in the high contrast limit equation (jS)) reduces to Q 



E(x) - Eb(x) ^ [ 



,^(xWO 
47r|x - x'l 



V^(x), (4) 



in which the potential term V(p, which depends on the 
detailed form of Kb, is explicitly curl free, and it follows 
from ^ that the divergence of the first term is of relative 
order kI/k^ It follows as well that ^ enforces the 
boundary condition Eint • n = 0{eb/e) — >■ 0, where fi is 
the unit surface normal. The background E;, is to be 
treated as a known input here. 

Remarkably, there is a straightforward procedure for 
solving Q that entirely avoids computing ip, or the non- 
inductive part of Eb. Thus, let {ZAr(x)}^^^ be a com- 
plete set of basis functions supported on Vc, and obeying 
V • Zn = 0, X e Vc, and n • Zjv = 0, x g dVc- One may 
then expand 



^«(x)2Ei„t(x) = ^CwZAr(x). 

N 



(5) 



The key observation is that the inner product Jy (fxZn ■ 
V/ = for any scalar function /. Therefore, inserting 
([5]) into (lU and taking the inner product on the left with 
one obtains the matrix equation 



{tiuo - c)c = Cb, 



(6) 



where the self-adjoint arrays O, C are defined by 
(fx- 



Omn 



3 Zm(x)*-Z^(x) 



47r /i(T(x) 



'^x / rf3^'ZM(x)*-Z^(x') 



d 



4tt\x — x' 



Cb,M = I d3xZM(x)* •Eb(x). 



(7) 



Given the solution of this equation for the amplitude ar- 
ray C,, the result of an external magnetic field or EMI 
measurement is also independent of ip. The curl in the 
magnetic field relation B = (ia;)^^V x E annihilates V(^, 
while an induced voltage measurement involves a line in- 
tegral around the closed receiver loop Cr, which is also 
insensitive to any gradient contribution. 

Now, in the time domain, one is primarily interested 
in freely decaying signals when S,E;,, hence C,b, vanish. 



One seeks solutions in the form of a mode expansion 
E = AnE*^"^ (x)e~'''"*, in which each (normalized) 
mode shape E^") is expanded in the form ([5]), and the 
decay rates A„ and corresponding basis function coef- 
ficients C*-"-* are solutions to the generalized eigenvalue 
equation. 



(8) 



Using mode orthogonality it can be shown that the exci- 
tation coefficients An are given by 



r(«) 



E(")* • d\, L 



in) 



dte^-'dtlT{t), 



(9) 

in which Ct is the closed transmitter loop, and Irit) is 
the transmitter pulse, terminated at t = 0. The mea- 
sured EMI voltage is also a multiexponential sum, with 
amplitudes given by similar line over C^: 



V, - A, 



(f E(") • dl. (10) 



For general targets, explicit forms for the Z^r may be 
hard to come by, and for general (t(x), the integrals in ([7]) 
must be done numerically. However, for uniform, ellip- 
soidal targets both of these problems are absent. First, 
let zl^^^'' be basis functions for the unit sphere. We will 



= 1,2, defined by ^ 

ZVJ^i^) = r'+2PX„„(0,^) 
,„^(x) = V X [(1 



use the forms zj^^^^, i 
.(1) 

zp) 



^)r'+2fX,„(0,0)] (11) 
where X;,„ are the vector spherical harmonics [l|, and 



p = 0,l,2,.... For an ellipsoid the forms 

ZAr(x) = y ]aaZ'-^^^\xi/ai,X2/a2,X3/a3)eo 
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(12) 



have the desired properties, where are unit vectors 
along the principal directions. 

The key property is that Z^^^' are polynomials (of de- 
gree l + 2p), and hence so are Z^r, and, as advertised, the 
Chandrasekhar approach Q allows Coulomb integrals of 
the form in ^ to be performed analytically. The inner 
product of this result with another basis function, as in 
([7]), is then trivial to compute. The method is elegant, 
but somewhat intricate, and so will only be outlined here. 
The basic result we use is @: 



</.(x) ^ Jd^ 



„,p[m(x')] 



dt 



VXi)-^[/i(t;x)] 



(13) 



where p is any ID function of the quantity /i(x) = 
J2a^a/0'a ^ [O7I], which traccs out the family of sim- 
ilar ellipsoids inside Vc, i^{p) = p{p!)dij! is the an- 
tiderivative of p, p{t,y.) = Ylia^tli'^t + t) € [0,1], and 
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FIG. 1: (Color online) Three lowest order vertically circu- 
lating mode shapes, and corresponding decay rates, for a 
10 X 10 X 40 cm radius aluminum prolate spheroid: {Ey,Ez) 
are plotted in the x = plane. 



t(x) = for X E Vc, and is otherwise the solution to 
1 = J2a=i ^a/ (.'^a + '''): ^ ^ ^c, which defines the family 
of confocal ellipsoids surrounding Vc- 

We apply this relation to compute potentials 0„ due 
to monomial densities of the form Pn{n) = (1 — 1^)" /n\, 
so that Vn(l) — i'nifJ') = Pn+i{p)- The density in (fT3|) 
is then a polynomial of degree 2n. To isolate a single 
monomial, one takes a combination of derivatives with 
respect to x and a~^: 



D 



21+m 



(x) 



^ ^2k+m 



(14) 



where 1 = {liTh^h) is a vector of nonnegative integers, 
m ~ (toi, TO2, 'tis) with each TOq = or 1, |1 + m| = 
Eailc+m^), di = UJ-'^y-d'^/dia-^y', the primed 
sum is restricted to Ikl < 111 -|- 1, and 



^(21+m), X _ 



(|1 



7raia2a3 ■r^' 
-|k| + l)!^ ^k+™+p(r) 



a 

where the primed sum is over all < Pa < la^ and the 



a^^ = 1 cm 
ci = 3.5x ib'S/m 
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FIG. 2: (Color online) Decay rate spectrum vs. aspect ratio 
a = az/uxy, with fixed conductivity and radius ttxy The 
first 125 decay rates (computed using 232 basis functions) 
are plotted for each 0.1 < a < 10. Blue dots at a = 1 
show exact analytic results for the sphere, with degeneracy 
effects from enhanced symmetry evident. Red dots at a = 4 
mark the three modes shown in Fig. [T] (three other nearby 
modes circulate around the cylinder axis, rather than along 
it). Green dot at the right is the analytic result for the lowest 
mode for an infinite cylinder, a — ^ oo. More slowly increasing 
branches at the left, a <^ 1, correspond to modes with current 
patterns circulating in the xy-plane. 



combinatorial factor is 

^(feM ^ (-1)"+' ( I \ r(|-m-p)r(fc + m + p + l) 
^ kl \pjT{^-Tn-l) r(fc + m+i) 



Finally, the elliptic-type integrals are defined by 

dt 



Mr) 



(16) 
(17) 



These obey the iterative relation dA^/da^ — —{k^ + 
^)Ak+e„i and evaluation may be reduced to that of 
A(a, r) = Aq. With the convention ai > a2 > 03 one 
one obtains A(a, r) = 2F{9, q)/ — a|, in which F is 
the elliptic integral [lH, sin2(6') = yj{al - a§)/(a? + r) 
and = {a\ — a2)/(«i ~ ^i)- Derivatives of F, and the 
auxiliary elliptic integral E{ip, k) [lH, may be expressed 
back in terms of E, F, allowing ([T7|) to be evaluated by 
iteration. In the case of spheroids, ai = a2 or a2 = 03, 
fully analytic expressions for E, F are available (ll| . 

Numerical implementation proceeds as follows: (a) Use 
the polynomial forms (fTTj) . (IT2]) and Coulomb integrals 
(fT4l) . (fTSjl (with r = 0) to assemble (truncated) arrays 
O, C. (b) Solve the eigenvalue equation ([8]) for the decay 
rates and (internal) mode shapes, (c) Use the known 
transmitter geometry and pulse waveform to compute the 
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FIG. 3: (Color online) Comparisons of data (solid lines; taken 
by the NRL TEMTADS platform [3) and theoretical pre- 
dictions (dashed lines) for a 5 x 5 x 10 cm radius aluminum 
prolate spheroid (top) and a 10 x 10 x 4 cm radius aluminum 
oblate spheroid (bottom), at various depths-to-center dc, with 
various axial tilt angles 6. The multipliers indicated in each 
legend entry reflect a ~ 10% variability in the transmitter cur- 
rent, and are applied to the data to optimize the fit. Straight 
dashed lines show the predicted early time divergence 

[^. This first principles agreement, over nearly two decades 
in both time and voltage, is remarkable. 

excitation coefHcients ([9]) . This involves an external field 
computation, the Coulomb integral of Eint in Q , which 
follows from ([H]), with r > 0. (d) Finally, use the 
known receiver geometry to compute the voltage (|10p . 

Figure [T] shows three of the computed mode shapes, 
obtained using the 232 basis functions ([TT|) correspond- 
ing to I < l + 2p < 7. Figure [2] shows the computed decay 
rate spectrum for a two decade range of spheroid aspect 
ratios. Figure[3]describes comparisons with experimental 



data from aluminum spheroids. The truncation loses ac- 
curacy at early time since more rapidly decaying modes 
with more complex geometry are not properly captured. 
However, the predictions merge smoothly with those of 
the complementary early time theory Q, which predicts 
a divergence. The combined result yields an excel- 
lent fit over the full dynamic range. With this level of 
agreement (afforded equally by the present theory and 
the remarkable hardware improvement the subtle 

changes in curve shape with depth and orientation can 
indeed be inverted for target geometry. This will be dis- 
cussed in detail elsewhere. 
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[2] Ground penetrating radar operating at GHz frequencies 
provides some imaging capability, but works well only 
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ducing finite ip even in the limit kj/k — ^ 0, plays no role 
in what follows, and so will not be discussed further. 
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these fail to diagonalize the system (|6]) or (|8} :4:j. Their 
desirable properties are more than compensated by the 
analytic simplicity of the spherical harmonics. 
[11] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, 
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